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ABSTRACT 


A method for determining whether a wavelet is 
minimum-delay is evaluated; the method is used to show 
that the impulse response of the seismic recording 
instruments and seismometers used at the University 
of Alberta in 1971 is minimum-delay. 

The transfer function of the seismic recording 
instruments and seismometers is used to generate the 
impulse response and the amplitude and phase spectra. 

A least-squares inverse to the impulse response is 
convolved with the seismic reflection data to produce 
instrument deconvolved data. 

The method of least-squares predictive deconvolu- 
tion is discussed. Single channel and multichannel 
predictive deconvolution methods are tested on some 
seismic reflection data and the results are compared. 
Both techniques seem about equally effective in 
resolving wave trains into approximations to reflection 


coefficient time series. 
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CHAPTER 1 


INTRODUCTION 


1.1 Crustal Refraction and Reflection Studies 


The crust, or outer layer of the earth, repre- 
sents less than two percent of the earth's volume, 
yet it dominates man's view of the planet. Of primary 
importance as the source of economic wealth in the form 
of petroleum and mineral deposits, it has long been 
studied in an effort to locate such riches. As near 
surface deposits are discovered and exhausted, searchers 
must probe more deeply into the crust, The details of 
its structure as deduced by all geological and geophys- 
ical methods become increasingly more important, In 
addition, man studies the crust in an effort to predict 
such catastrophes as volcanoes and earthquakes and to 
minimize the damage they cause, Although such phenomena 
originate below the crust, we can only observe them 
through the crust, and we must know how the crust affects 
such indirect measurements. 

“the first indirect measurements of crustal struc- 
ture were made using near earthquake observations. In 
1910, a velocity discontinuity was discovered by 
Mohorovicié,and it has been defined as the base of the 


crust. It has been debated whether the discontinuity 
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represents a phase or a chemical transition, In 1925, 
a shallower velocity discontinuity was discovered by 
Conrad. It has been regarded as a non-universal tran- 
Sicion from granitic*to basic material. 

The first field tests of the reflection seismic 
method applied to exploration were conducted in Oklahoma 
Mito 2 le svochricver,) 1952). the first recording of deep 
crustal reflections was made in Montana (Junger, 1951), 
where reflections were recorded in the time interval 
7.0 to 8.5 seconds. Reports of deep crustal seismic 
reflections have been documented (Steinhart and Meyer, 
1961). A good review of deep seismic sounding in 
northern Europe has been presented (Vogel, 1971). The 
American Geophysical Union series of geophysical mono- 
graphs has periodically reviewed seismology (Steinhart 
and Smith, 1966; Knopoff, Drake, and Hart, 1968; Hart, 
1969; Heacock, 1971). 

A review of studies of crustal structure in 
Germany has been published (Dohr and Fuchs, 1967). 

A statistical study of deep reflections was done; 
histograms were plotted with reflection time versus 
number of reflections for a small region. It was found 
that reflections existed on the records in a statistical 


sense. The deep reflections were correlated only over 


short distances. 
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It has been suggested that the Mohorovi¢ié dis- 
continuity is stepwise, interrupted by a rhythmically 
arranged series of partial melts of lower velocity 
(Meisner, 1967). It has been shown that a simple layered 
model of the reflecting horizons is not consistent with 
the observed large amplitude and low cut-off frequency 
of deep reflections (Fuchs, 1969), but that a transition 
zone with a series of velocity reversals provides a 
satisfactory model, 

The first seismic crustal study in Alberta was 
undertaken by a group of fifteen seismic parties 
(Richards and Walker, 1959). Measurements were recorded 
along a 130 kilometer reversed refraction profile, paral- 
lel to and 100 kilometers east of the frontal thrust of 
the Rocky Mountains. The depth to the Mohoroviéié dis- 
continuity was found to be 43 kilometers, with an upper 
mantle velocity of 8.2 kilometers per second. An inter- 
mediate discontinuity was found at a depth of 29 kilome- 


ters with a velocity of 7.2 kilometers per second below. 


1.2 Crustal Studies at the University of Alberta 
Explosion seismology was first used at the 

University of Alberta when a 250 kilometer east-west 

reversed refraction profile was shot in southern 


Alberta (Cumming, Garland, and Vozoff, 2962) ae 
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crustal model was obtained which consisted of three 
layers between the basement and upper mantle. The 
refraction work was extended along two east-west 
reversed refraction profiles for 500 kilometers 
(Cumming and Kanasewich, 1966). In addition, a 
refraction profile extending west into the Rocky 
Mountains was shot. 

The first deep reflection work at the University 
of Alberta was a study of near-vertical incidence 
seismic reflections at Lomond, Alberta (Kanasewich and 
Cumming, 1965). Good reflections were recorded at a 
time of 11.4 seconds along an expanding spread. The 
use of three-component geophone arrays showed that the 
event was near-vertical incidence. A 7 versus — 
analysis gave an average velocity from the top of the 
Mississippian to this reflector of 6.37 kilometers per 
second and a depth to the reflector of 34 kilometers. 
The discontinuity has been referred to locally as the 
Riel discontinuity to avoid implying that the Conrad 
is universal. 

Seismic reflections were studied over four pro- 
files for a total of 90 kilometers (Clowes, Kanasewich, 
and Cumming, 1968). Autopower spectral calculations 
showed that the energy of the reflected wavelets was 


concentrated in the five to fifteen Hertz range. Along 
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one profile, the Riel reflection was correlated over 
nearly .25 kilometers... Analysis,of this data.gave a 
revised average velocity of 6.2 kilometers per second 
to a depth of 34 kilometers. 

The attenuation properties of the crust and 
the nature of the reflecting horizons were. studied 
(Clowes and Kanasewich,1970). An attenuation model 
was obtained from a comparison of the data with syn- 
thetic seismograms which included attenuation as a 
function of frequency and depth. Autopower spectral 
analysis was used to examine the nature of the reflect- 
ing zones. It was concluded that models with first 
order discontinuities, linear gradients, and step-like 
velocity increases did not match the spectral charac- 
teristics of the reflections. A transition zone model 
consisting of sills of alternating high and low velocity 
material was in agreement with the observed data. The 
sills of the model were less than 0.2 kilometers thick 
extending over a zone less than 1 kilometer thick. 

Studies of crustal. structure.in; central Alberta 
were undertaken using P-codas from earthquake seismo- 
grams, (Ellis and, Basham,,1968)...Vertical. to, radial 
spectral ratios were compared with theoretical ratios 
calculated using the Thomson-Haskell matrix formulation, 


With the aid of synthetic seismograms, an improved model 
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was obtained (Somerville and Ellis, 1972). Another 
crustal study in the same area was done using the P- 
coda spectral ratio method (Sprenke, 1972). 

Synthetic seismogram calculations using asymp- 
totic ray theory for plane-layered, homogeneous elastic 
media have been described (Hron and Kanasewich, 1971). 
Criteria for the selection of phases which make signi- 
ficant contributions to synthetic seismograms have been 
given (Hron, 1971). The synthetic seismogram calcula- 
tions have been generalized to permit curvilinear inter- 
faces in the media (Hron, personal communication). 

A 760 kilometer refraction profile from Greenbush 
Lake, British Columbia to Swift Current, Saskatchewan 
was interpreted (Chandra and Cumming, 1972). A unified 
interpretation was found from an examination of previous 
refraction and reflection data and from gravity and 
magnetic surveys. 

Investigation of an additional 40 kilometers of 
reflection data in southern Alberta has been undertaken 
(Cumming and Chandra, personal communication). The 
studies have used autopower spectral analysis and homo- 
morphic deconvolution of the major reflecting zones of 


the deep crust in order to further define the details 


of the reflecting interfaces. 
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A seismic reflection crustal model has been 
obtained near Edmonton, Alberta (Ganley, 1973) 
Velocity analysis of a profile 13 kilometers long 
has indicated that fifteen degree southeasterly dips 
occur within the crust, with a total crustal thickness 
of 35.5 kilometers. 

Homomorphic deconvolution has been used in 
analysis of earthquake seismograms (Alpaslan, personal 
communication). 

Much of this work has indicated that deep crus- 
tal reflecting and refracting transitions are zones of 
alternating high and low velocities of less than 1 
kilometer in extent, and the resolution of normal 
seismic records is insufficient to permit examination 
of the detail of these zones due to the overlap of 
pulses. Therefore, it was thought that it would be 
worthwhile to apply the techniques of deconvolution 
to reflection records from the deep crust in an effort 


to obtain more information about the structure of the 


transition zones. 
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CHAPTER 2 


CONVOLUTION AND MINIMUM-DELAY 
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Seismic reflection data is recorded in either 
analog or digital form. An analog trace is a con- 
tinuous physical variable, while a digital trace is 
a discrete physical variable. A digital trace may be 
represented by a series of complex numbers which are 
samples or measurements of the physical variable at a 
constant time interval At. Thus, the trace X may be 
represented by the series KX yrXprXoree- corresponding 
fo -Sampaes Otters at times 0, At) 2A, Ge.u0 co 

The mean of X is defined to be 


1 N 
lim N ) x ‘ e oD, x i. 
N00 j= 


The autocorrelation of X is defined to be 


where the horizontal bar means the complex conjugate. 
The trace E may be defined as white noise (Robinson, 


1964). Then E must satisfy the following: 


geddie ni bebvooss al sdsh!noltoeltes olmesse | 
“109 6 ek eDBIt poise oA mro2 istiptbh 10 poisnas 
al sostd {Btip2b 5 aad .eldsixsy fsoteydg evounts 
4d ysm eDat? ted cor) A .oldsizev Isoieydq storoesb 5 
sis foidw aredmun. ¥yelqmoo Io esiise s yd bsgnessigez 
5 #6 Sidstisv Ispieydq odd to etnomeivesom 10 e9lqnse 
ad yam X sosit oft ,eud? .3A Isvasdai ombt gasteaol — 
pribnoqesix0d capes peig® 28t382 sit yd bednezssges — 
+ ewe (HAS. 2a) 90 BOM te K 39 eolguse oF 

od of beaiieb ai X Fo asst edD | 

wa < 


7 MM 
ts : x tT = mil 


& 
Ca 


we 
Pe ene eee 
ae 


1. E must have zero mean; 

2. the zero lag of the autocorrelation of E, 
Lor must be finite; and 

3. all other lags of the autocorrelation of E 
must be zero. 


The seismic trace may be defined to be the con- 


volution of white noise Cpr epregress with a character- 
istic wavelet By rb bor-- rb, by the convolution 
equation 
co 
sae ) Dano Bo 
s=0 


This equation may be written symbolically as 


X = BYE. 2.4 


The object of deconvolution is to obtain E given XxX, 
since E represents the reflection coefficients of the 
earth beneath the seismic instruments - the reflection 
coefficients give details about the structure of the 
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A two-length wavelet by by has a reverse Shpeaiee 
The minimum-delay wavelet of the pair is the one in 
which the first sample is largest in magnitude; the 
other wavelet is maximum-delay. The convolution of 


minimum-delay wavelets gives a minimum-delay wavelet; 


the convolution of maximum-delay wavelets gives a 
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maximum-delay wavelet; and the convolution of one or 
more minimum-delay wavelets with one or more maximun- 


delay wavelets gives a mixed-delay wavelet. 


252m Zexos of:a Polynomial 


The Z-transform of the wavelet Bg by rbor-- eb, 


1s"ar polynomial“G@n*'z"Got order n> 


2 n 
bo + b,2 + b,z ry es bi2 ; aS 


If we factor this polynomial, we obtain it in the form 
b. (2-21) (2-25) wees (z-z_) > 2.16 


The zeros of this polynomial are at z = ZyprZoreeer Ze 


A two-length wavelet Bo rb, has a Z-transform bot biz 
which can be factored to give b, [2-(-by/b,) 1. The 
zero of this polynomial is at Zi = “by /b, - If the 
wavelet was minimum-delay, [by | 2a (Dy [ye SON Mee eee 
If the wavelet was maximum-delay, [Bo | < Ib, |, so 
|2,1< 1. The boundary |z| = 1 separates minimum-delay 
zeros from maximum-delay zeros; this boundary is the 
unit Circle in the complex z-planes If we plot the 
zeros of a polynomial in the complex z-plane, we may 
determine whether the wavelet is minimum-delay or not. 


Thus one method of determining whether a wave- 


let is minimum-delay or not is the factoring of the 
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polynomial obtained by Z-transforming the wavelet. If 
the wavelet is longer than a few samples, a general 
method of factorization is needed. Subroutines have 
been written to factor a polynomial; subroutine DPROD 
from IBM's Scientific Subroutine Package was found to 
be the most satisfactory. DPRQD uses double-precision 
arrays, and it does not limit the order of the input 
polynomial. It provides an error return code which 
indicates any existing error condition. The subroutine 
takes considerably less than one minute of Central 
Processing Unit time on the IBM 360/67 at the University 
of Alberta Computing Center to factor a.polynomial,. of 
order forty. 

To test DPRQD a trial wavelet was chosen, the 
expansion of (a1 .95)>7 (z-1.1)?,which has been plotted 
im Figure 251. They resulting polynomial’ of order ftourteen 
has a zero of order twelve at z = -1.75 and a zero of 
@rder two at z= 1.14 The zeros calculated by DPROD 
have been plotted in the upper half of Figure 2.2. 

Note that only the upper half of the complex z-plane 
is plotted, since zeros occur singly on the real axis 
Or as complex conjugate pairs symmetrical about the real 
axis. While the zeros show a large scatter, roughly on 


a circle centered at z = -2.0 with radius 1.0, it may 


u | ae 


32 .delovew itiaeiieninns eileen” | 
feionse © Zeke WeE  nhilt aRRHOR'S Sotoent SS 7 
svsd esnissosdue .bebser ef nottssixédoni. to bodsem 
doada sattvercve 4 -Isimonylog 6 - ‘totos? et asttixuw asd 
+ bavot esw sps2zost oniisordue olgeseetee) a'Mai moxt 
iit eoey GOAGG .yrotosteitse daom orit od 
tuqni sig to +9510 Se dimiD-som, 2290p, +h, bus ,eysu7s 
dotfw shoo mivisr totie m6 eebivorg at . fsimonylog 
enitvoxudue oAT «nots Ebnoo yoris patseinxne ys estanibai 
fe13n99 to satunift Sa0 nes 2zeol yidsxobi anos eodst 
yiiexevinu oft 36 Ta\0ae MAL ois no omis tin pateesoord 
to Isimonylog 6 tos5st 95 ‘164099 prisuqmod ssrsdté to 
~Yt102 tebee 

ei3 \aeeoio esw toiovew {pixtt 6 aguda te83 of 
bestola aesd sed dsidyss (1.f-s) ent 32 f+) to roisnsgKe 
destuo? aebio to Isimonylog paidivesr sat -t. § suvpit ak i 
to oss 6 Bas c¥.f[- = = 36 ovisw2 +sbi0 To often 6 esd 
dqQada yd bedslvolso eouss Si? .1,f = s 48 owt tebz0 7 
~S,S saveit to tisd xeqqu aft at bedtolq reed oved of 
eusiges xoiqmoo end to Fisd zsqqu odd yino tens eae ai 
aixs Less odd no. yipate ivo90 20198 sonte \bestolg a 
teow est) suede Seoztemmye cxisq easputnon xelqnos as 30 
ao Yideuos «xettso2 optst & wore goxes oft of kaw wales 
= +2 0.0 eutbex datw 0 .S- = 5 $5 betetmes eforto 6 


FZ 


0.60 1.00 


0.20 


6 
SAMPLE 


AMPLITUDE 
~0.20 


-0.60 


-1.00 


Figure 2.1. The trial wavelet. 
Amplitude units are plotted versus sample units; 
the amplitude has been normalized to 1.0. The 
trial wavelet is minimum-delay. 


153 


Figure 2.2. Zeros of the polynomials. 
Upper diagram: zeros of the trial wavelet. 
Lower diagram: zeros of the impulse response of the 
seismic recording instruments and seismometers used in 
1971. In each case, only the upper half z-plane is plotted. 
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be noted that all the zeros are outside the unit circle 
as required. 

DPRQD was also tested on a longer wavelet. The 
impulse response of the seismic recording instruments 
and seismometers used at the University of Alberta in 
1971 was generated. The seismic recording instruments 
and seismometers are described in Chapter Three. The 
impulse response was sampled at an interval of At = 
0.0056 seconds, and it has been plotted in Figure 2.3. 
It was truncated after twenty-four samples and Z-trans- 
formed; the polynomial of twenty-third order was input 
to DPRQD. The zeros of this polynomial have been plotted 
in the lower half of Figure 2.2. It may be seen that 
twenty zeros fall outside the unit circle, roughly on 
a circle centered at z = 0.0 with radius 1.2. Two other 
zeros fall on the unit circle where’ it crosses the posi- 
tive real axis. The final zero falls on the negative 
real axis inside the unit circle. Since the mean of the 
samples was not removed, it may be suggested that remov- 
ing the mean would put this zero on or outside the unit 
Circle (O'Brien, 1969). it may be concluded that the 
impulse response is minimum-delay. This agrees with the 
conclusion drawn from examining the transfer function 
of the instruments and seismometers in the complex s- 
plane. The transfer function of the seismic recording 


instruments and seismometers is described in the Appendix. 
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Frqure. 2,3. The impulse response of thevseismic 
recording instruments and seismometers used at 
the University of Alberta in 1971. The sample 
interval At is 0.0056 seconds. 
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DPRQD has an error return code which indicates 
whether or not the expansion of the calculated zeros 
is identical to the input polynomial to six figures of 
accuracy. Although the calculated zeros of the trial 
wavelet satisfy this condition, they are not the same 
as the zeros which were originally used to generate the 
trial wavelet. Shorter wavelets, which were expressed 
exactly by six figures of accuracy, were input to DPRQD, 
and their zeros were recovered exactly. This indicates 
that the zeros calculated for the trial wavelet are 
correct, and that the reason they show scatter is that 
the trial wavelet must agree with the expansion of the 
calculated zeros to only six figures of accuracy. In 
conclusion, although the zeros found by this procedure 
show scatter, they appear to be sufficiently well resolved 


to indicate whether or not the polynomial is minimum-delay. 
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CHAPTER 3 


INSTRUMENT DECONVOLUTION 


3.1 Seismic Recording Instruments and Seismometers 


In 1969, the University of Alberta seismic 
reflection equipment was deployed in two trucks; in 
each truck, one Texas Instruments VLF-2 amplifier 
system recorded six channels of data on a Precision 
Instruments seven track FM analog tape recorder. The 
s-plane representation of the seismic recording 
instruments is given in the Appendix. The impulse 
response of the seismic recording instruments is 
plotted in Figure 3.1, and the amplitude and phase 
spectra are plotted in the Appendix. 

In 1970, the University of Alberta seismic 
reflection equipment was deployed in one truck, in 
which one Texas Instruments VLF-2 amplifier system 
recorded twelve channels of data on a Peripheral 
Equipment Corporation model 7830-9 digital tape 
recorder. The s-plane representation of the seismic 
recording instruments is given in the Appendix. The 
impulse response of the seismic recording instruments 
is plotted in Figure 3.2, and the amplitude and phase 


spectra are plotted in the Appendix. 
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Figure 3.1. The impulse response of the seismic recording 
instruments used at the University of Alberta in 1969. 
The sample interval At is 0.0056 seconds. 
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PPGULe #72) sl he Impulse” response OL Che seismic recording 
instruments used at the University of Alberta in 1970. 
The sample interval At is 0.0056 seconds. 


bn) 


- 


7 


iv 
- . 


he 
= 


- 


oa 3 
—— 


= 


_ 


20 


In 1971, the University of Alberta seismic 
reflection equipment was deployed in one truck, in 
which one amplifier system (Alsopp, Burke, and 
Cumming, 1972) recorded twelve channels of data on 
a Peripheral Equipment Corporation model 7830-9 
digital tape recorder. Each channel represented the 
output from a tapered array of sixteen Electro-Tech 
EVS-4 seismometers (Clowes, 1969). Two additional 
channels were recorded; one channel was shot instant 
from a radio receiver, and the other channel was 
absolute time from WWVB. These fourteen channels were 
recorded on a nine track synchronous tape recorder at 
a density ‘of 800 bytes per inch with a tape velocity 
of 6.25 inches per second, resulting in a data transfer 
rate of 5000 bytes per second. 

The signal was converted from analog to digital 
form by’ a-converter*with a resolution of thirteén bits 
plus*Sign providing=a dynamie*range*or eighty-four 
decibels. Two eight-bit bytes were required to store 
each sample, so the sample rate was 2500 samples per 
second. Each channel was therefore sampled once every 
14/2500 or 0.0056 seconds, and the Nyquist frequency 
was slightly above eighty-nine Hertz. The length of 


the records written was 4096 samples for a time inter- 


val of 22.94 seconds. 
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The multiplexer, analog to digital converter, 
and five thousand Hertz crystal clock plus associated 
interfaces and controls were contained in a modified 
model 120 Data Acquisition System manufactured by 
Datum Inc. Controls allowed the selection of from 
one to twenty channels to be sampled sequentially 
and the selection of record length up to a maximum 
of 4096 samples per channel. 

The amplifier circuit had.a flat response from 
one-tenth to fifty Hertz, and it consisted of an 
input transformer of eighteen decibels gain followed 
by a preamplifier of forty decibels gain connected to 
the amplifier by a two pole high pass Bessel filter 
which could be bypassed. The amplifier had switchable 
gain, and its output passed through a four pole low 
passevessel [filter with aesiiity Hertz core. 

The s-plane representation of the seismic 
recording instruments and seismometers used in 1971 
is given in the Appendix. The impulse response of 
the seismic recording instruments and seismometers 
is plotted in Figure 2.3; the impulse response of 
the seismic recording instruments is plotted in Figure 
3.3. The amplitude and phase spectra of “the seismic 
recording instruments and of the seismic recording 


instruments and seismometers are plotted in the Appendix. 
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Figure 3.3. The impulse response of the seismic recording 
instruments used at the University of Alberta in 1971. 
The sample interval At is 0.0056 seconds. 
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3.4 instrument Deconvolution of the Data 


The transfer function in the s-plane of the 
seismic recording instruments and seismometers was 
obtained in the first step of See pane ae wemeoahes 
tion. An expression for each stage of the seismic 
recording instruments was calculated using Laplace 
transform theory and laboratory testing; the general 
s-plane expression for the seismometers is derived 
in the Appendix. 

A computer program was written by Cumming which 
calculated the amplitude and phase spectra, as well 
as impulse response, for a given transfer function. 
The program used a basic Fourier transform subroutine, 
FASTFT, which made use of subroutines COSP, COSTAB, 
SINTAB, and SPLIT (Robinson, 1966). The amplitude 
and phase spectra were calculated over a frequency 
range from one-hundredth to one hundred Hertz; the 
impulse responses were sampled at an interval of 
0.0056 seconds for 1.25°seconds; and’ all plots*were 
done using a model 770/663 CALCOMP Plotter at the 
University of Alberta Computing Center. 

The impulse response of the seismic recording 
instruments and seismometers used at the University 
of Alberta in 1971 was generated and truncated after 


forty-three samples; it was then input to a program 
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which used subroutine SPIKE, which made use of sub- 
routines IMPULS, SHAPE, CROSS, EUREKA >> DOT, sFOLD; 
ZERO, and MINSN (Robinson, 1967). SPIKE was used 

to compute a spiking filter which could be convolved 
with the truncated impulse response to give an 
approximation to a delta function. The spiking 
filter was then convolved with the data to give 

data which had been instrument deconvolved. 

Let Do rby rbo1--- rb, be the truncated impulse 
response, and let dp rdyrdor--- sd be the desired 
output, a delta function, which has all its coefficients 
equal-=Eowzeronexcept atsangle’ coebiicient’ which is equal 
tovones) Ihe purpose of “SPIKE is to find taspiking 


filter f 6 2et such that when we convolve the 


oe aes 

truncated impulse response with the spiking filter, we 
obtain an output wavelet CorCy Corser Cn which is as 
close as possible to the desired output in a least- 
squares sense. 


Tiev spiking filter “is calculated ‘so (tliat the 


sum of squares given by 


is a minimum. Substituting for Cy, we have 
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The value of I is a minimum if its partial derivatives 
with respect to each of the coefficients Eqrtyrtareees 


En are vequal "to *zero.eiThusiwe have 
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The autocorrelation of the truncated impulse response 


is defined as, using equation 2.2, 
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cross correlation between the desired output and the 


truncated impulse response is defined as 
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This set of m+l simultaneous linear equations in the 
unknown coefficients a a is the discrete 
time analog of the Wiener-Hopf integral equation 
(Robinson, 1967). The equations are called the normal 
equaeigns Loretnieeti tter: 

The spiking filter was arbitrarily chosen to be 
of length forty-three, since the least-squares error 
decreased monotonically as the length of the spiking 
filter increased. SPIKE was used to calculate the 
léast-squaresSseYTroretom Sacheof=thespocssible positions 
of the one in the delta function; the error was a 
minimum when the one was in the twenty-ninth position. 

The spiking filter was convolved with the data 
obtained at the University of Alberta in 1971. The 
OnLogagnededata, sf iiitered, bye daeiuve@ tos Liaely Hertz zero 
phase shift eight pole Butterworth bandpass filter 
(Al pds ane 6.3) Fees oOUreda iahguLege. 4. Fihe data, 
convolved with the spiking filter and filtered with the 
same bandpass filter as before, is plotted in Figure 
3, Se Althougmethestwoseplotsfareavery samilar,; there 
is a definite sharpening of several three sided pulses 
into two sided pulses. The resulting two sided pulses 
resemble the half cycle sine wave pulses obtained for 
the Pierre Shale (McDonal et al., 1958). The instru- 


ment deconvolution is thus moderately successful. 
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Figure 3.4. Eleven channels of data from the first 
shot on September 2, 1971 before instrument 
deconvolution. 
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Eleven channels of data from the first 


shot on September 2, 1971 after instrument 
deconvolution. 
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CHAPTER 4 


PREDICTIVE DECONVOLUTION 


4.1 The Predictive Decomposition Theorem 


The following discussion of the Predictive Decom- 
position Theorem and of the method of predictive decon- 
volution is based largely on three references (Robinson, 
1954; Robinson, 1967; and Peacock and Treitel, 1969). 

Any observational time series x, (-~» < t < o) 
may be considered as a realization of a random process, 
which is a mathematical abstraction defined with respect 
to a probability field. A time series is said to be 
stationary if the probabilities involved in the random 
ProcessedOsmoOl nave a Speci iiceora gate ti atame. 

For any random process, averages may be taken with 
respect to the ensemble of realizations X, for a fixed 
time’ t; such averages are called ensemble averages. The 
ensemble average or mean value of x, is denoted by E[x,], 
and the variance is denoted by EL (x, - E(x,1)71; both 
the mean and the variance of a stationary random process 


are independent of time. 


The autocorrelation coefficients 
= b Aol 
: (t) B(x. 1 4%,] 


are independent of s, and they are an even LUNCCLON LOL 
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the time t, that is 


Dy (t) = UB peteh=), 5 gets) 


Another type of average is the time average, 
taken with respect to all values of time t for a fixed 
realization xX (-ocait <.%) Of the ‘random process. 
Stationary process us called an ergodic process if the 
ensemble averages and the time averages are equal with 
probability one. As a result, the autocorrelation of 
an ergodic process may be expressed as a time average, 
as in earlier chapters. 

The autocorrelation is a non-negative definite 


function; that is, equation 4.2 holds, and 


iH] 
bt 
Nh 
= 
bo 
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n n 
a) a Peed aay Se n 


for every real set of Ay 1An Agree ery (Robinson, 1954). 
This property of the autocorrelation is equivalent to 


its representation by the Fourier transform 


1 
(t) = = | GOS Wee autOr 4.4 
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a representation known as the Wiener-Khintchine theorem. 
The spectral distribution function A(w) is a real mono- 
tonically increasing function of w (0 s<w<« T) with 


NCOs OVand A(T), =sre, Lhe inversion formula is 
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BE pg | Pxx (+) | is convergent, then the spectral 
distribution function will be absolutely continuous, 
with the continuous derivative 


aA (w) 


O(w) = an 


where @(w) is the power spectrum, an even, non-negative 
function. In terms of the power spectrum, the auto- 


correlation is 


T 
b(t) = oF | exp (it) o(w) dw. 4.7 
T 


A time series with an absolutely continuous spec- 
tral distribution function is generated by a process of 
convolution. For a fixed realization of white noise 


, the corresponding fixed realization 


Fe eee 
of a process of convolution is See See ee where 
= -~< t< 4.8 
x ) Cee we for all (ce co 


It may be assumed that the e, are mutually uncorrelated, 
and that the mean and the variance of the e, are zero 


and one, respectively; then 


Ele.e,] = Ele.] Efe, ] =’ OP ror St. 4.9 
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The autocorrelation coefficients of a process of convolu- 


tion x, May be written 


(t) =E[ 
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By using the relationship 
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the autocorrelation may be written 
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j=-© k=-0 
= 
4.11 
WhLenh Ls in the torm of CGquation 4’. 7 where 
be an ae ; 
D (Gi= hCC) |=" ol Cc, exp(iwj) ) cy, exp(-iwk), 
j=z-© =—0 
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C(w) being the transfer function of the operator c,. 


Thus, the spectral distribution function is absolutely 
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continuous. Conversely, any process with an absolutely 
continuous spectral distribution function is a process 
Ofrconvolucion'’: 

It may be shown that a process with a power spec- 
trum ®(wW) may be made physically realizable, stable, and 
minimum-delay by using the Kolmogorov factorization of 
the spectrum or the Fejer factorization method (Robinson, 


1954). That is, the power spectrum may be factored 


6(w) = |B(w)|* = | J b, exp(-iwt) |? , Peat 
t=0 
where b, = 0) OG st es. 0 5 bo a ee 
y be oe; y |b, | <teo) (Tand 
t=0 t=0 


It may be seen that equation 4.13 may be used in place 
GOpmequation 4,12%+ i other words, the operator Co rly 1Coe 
228 May be replaced by the operator Dy rby bore Ay alg hults) 


the process of convolution is given by 


replacing equation 4.8 with its minimum-delay equivalent. 
In this equation, the Ce, (-o < t < ~) represent a real- 


ization of a mutually uncorrelated process with zero mean 
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andsunit wanilance»« thats white. noise. 

The Predictive Decomposition Theorem may now be 
stated. Given a non-deterministic stationary process 
x, (Sete <>) with -discretertimesiparameter t, suppose 
x, has an absolutely continuous spectral distribution 


ie: 


LUnCtLON. slhen there exists agdecomposttion 


AS boey os bie, 4 < b,e,_>5 UO 9 4.16 
where the Ce, have zero mean and unit variance and are 
mutually uncorrelated, and the operator Do rby bor... is 
minimum-delay with bo > 0 and 2 =. bt a bs Waa iar eal se Core 


This may be described as the minimum-delay specialization 

of the Wold Decomposition Theorem (Robinson, 1954). 
Because the operator Do rby bore. is minimum-delay, 

there exists an inverse operator AQ ray ragrees which is 


realizable and minimum-delay. Thus we have the inverse 


predictive decomposition of the time series xy given by 


GQ, = 4X, + a,Xp,_1 + 42%t_2 Tee es any 


4.2 Single Channel Predictive Deconvolution 

The model’ required for the application of predic- 
tive deconvolution is based on the hypothesis that the 
strengths and arrival times of the signals ina seismic 


trace can be represented as white noise, and on the 
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hypothesis that the characteristic wavelet associated 


with each signal is minimum-delay. 


nN 


The: prediction of Xp ty by Xi, is 


X a Pats®t-s = 2 Des ey ee ANE 


using equations 4.16 and 4.17. Now let r = stn; then 
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Now, define the filter operator a a where 
OGY are a) a 4.20 
520 ots ©r-s 
so that equation 4.19 becomes 
eta 7 FR (a)x, =f) la)x, +f, (a)x,_, +f, (a)x, 5+... 
4.21 
Thus, oe is expressed in terms of a linear combination 
of the present x, and the past X,_j/X,_5/X%,_37-++ + The 


operator fof, ,f5,--- is realizable since f(a) is zero 
for wr less than zero (Robinson, 1954). 
In practical applications, only a finite number 
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Nh EX ti fix] nh PoXp_9 +... +f x p 
4.22 
Using the method of least-squares, the value of 
nS y ( s ; 
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The value of I is a minimum if its partial derivatives 
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Theretore, the Lilter operaror Eqrfj fare - erin is given 


as the solution of these normal equations. The solution 
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is found by a recursive algorithm; if the filter operator 
is of length m, machine time is proportional to m*, while 
computer storage is proportional to m (Robinson, 1967). 


In matrix form, equation 4.26 is 
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The autocorrelation matrix on the left is symmetric, has 
its elements equal along each diagonal, has non-negative 
eigenvalues, and has a non-negative determinant. 


mie prediction error time yseries is, delined sas 
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From equation 4.24, the least-squares error is 
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Using equation 4.26 to evaluate the term in brackets, 


the least-squares error is 
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Equation 4.27, in matrix form, may be rewritten 


as a system of m+l simultaneous linear equations, 
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This system of equations may be augmented in such a way 
that the prediction operator is converted to its corre- 
sponding prediction error operator (Peacock and Treitel, 


1969). The matrix equation corresponding to the augmented 
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Thus, the prediction operator of equation 4.27 has been 


converted to its corresponding prediction error operator. 


For the case a 
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Equation 3.6 may be written in matrix form, Se een 
Go a Sali 0, G5 = 0% 50> OF me oss ty See F 0, recall- 
ing that the autocorrelation is an even function, and 


changing sm wtosm+1, 


ro ry r. 5 ee Ctl fo al: 
Ete oetned. xn fy 0 
ry ry ry eng | f. = 0 3 AG 37 
Tatl Tm tm-1 hi 0 Eel 0 


BUGetias equation is the one for the calculation, otetie 
least-squares inverse filter which compresses an unknown 
Signal into a unit spike at zero delay (Peacock and 
Preoueely, L969)°. Comparing equation 4.597 to equation 4.35, 
it may be seen that the m+2 length prediction error oper- 
ator with unit prediction distance is identical to the 
ZeLoude lay, least-squares inverse £1 1 Cer of sleng thi mi27 
except for a scale factor Bo: 

Similarly, equation 4.33 may be compared to 
equation 3.6; then hes Bor chee Bae bard Boe Seeley 
ce aT See ui oa Osta vey Sp Ree O See LNe SSOae LoneoOL 
equation 3.6, changing m to atm, is a least-squares 


inverse tiiter of wength otmtljechessolution (on equation 
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4,33 is the prediction error operator of length g+mtl 
for, the prediction distance of gq; therefore, these two 
filters are tdentical. 

Using equation 4.10, the autocorrelation of the 
input may be considered to be the autocorrelation of 
Liem NatiGCleiLouLlc, WAVGlCt. = bULCH LIM Glng tne: positive 
lag nalt of the’ autocorrelation to length g+mt+l implies 
that the characteristic wavelet is of length g+mtl. 

Since 5. = 0 for a < j < atm, and since a5 is 
the cross-correlation between the desired output and 
the input, the cross-correlation is zero for lags greater 
than ag-l. Therefore, it may be concluded that the desired 
output Is Gof length oo.” Since the desired output is of 
tencthieg, bts aucOCORLELALLON JS Zerorat lags. om, atl, 
note eee LUC PLCalCELlOl GriCneti i lel wilt. Modi tyethe 
Mipiesinestched way thd ue tle sal LOCOLre lation ole tne Out 
DUCEWLUeLena to ve 2zGLo dU slags q, otl, of2;7.... , Otm. 

THe predicuLoneerror Operacon Contracts Che char 
acteristic wavelet of length qtmtl into a wavelet of 
length q. Least-squares inverse filtering contracts a 
characteristic wavelet of length qgtm+l into a wavelet of 
length 1. Thus, it may be concluded that predictive 
deconvolution is a generalized method of least-squares 
inverse filtering in which the degree of wavelet contrac— 


tion may be controlled. 
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4.3 Test of Reverberation Removal 


Removal of predictable energy from data is 
possible if the method of predictive deconvolution 
is used; the energy removed will have periods between 
a and atm time intervals, and it may be classed as 
either short period or long period energy. The first 
step is to plot the positive lag half of the auto- 
correlation. Short period reverberations will appear 
as decaying waveforms on the autocorrelation plot, 
while long period reverberations will appear as distinct 
waveforms separated by quiet intervals. 

Removing short period reverberations is done 
differently depending on whether or not the decaying 
waveform on the autocorrelation plot is highly regular. 
If the waveform is highly regular, a is set at the 
desired prediction distance while atm+l is set to span 
Q plus one -completewcycle. sit theiwaveronm is Prreguiar, 
a is set at the desired prediction distance while atm+l 
is set to span a plus the remaining significant interval. 
In either case, a should be set to span the interval 
between zero lag and the second zero crossing on the 
autocorrelation plow (Peacock and Ereiltel = 1969)". 

Removing long period reverberations is done 
differently depending on whether or not the ringing 
TWsuotearlirst Orders Navurern) lim the: ringing. losor 


a first order nature, « is set to’ span the interval 
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between zero lag and the onset of the first multiple on 
the autocorrelation plot, while atm+l spans a plus the 
finest muLtpe. If the ringing tsvor ashigher order 
nature, a iS again set to span the interval between zero 
lag@and the onset of the first multiple on the auto-— 
COLeLatLion plot, wille dtm-lespans co, plus the first 

two or more multiples. In either case, repetitions of 
the waveform centered at multiples of a+ (m+tl1)/2 will 

be removed (Peacock and Treitel, 1969). 

In order to test the effectiveness of the method 
of predictive deconvolution, a trial trace was created. 
This trace used the trial wavelet of Chapter Two; the 
trial wavelet was convolved with the trace e,=1, e 


0 al 
Seated AU IRA UAnn C59= 0, C39= sky Aclelep aep cell’ (esaelersy: alte (edWene efarel 


= 0, 


on whe upper left of Figures 4.1,$4.2, and 4.3. The 
autocorrelation plot of the trial trace is plotted on 
Chemuipper srigit ot FPagires 4.1) -ang. 4.3. 

A main program was written to do single channel 
predictive deconvolution, using subroutine EUREKA 
(Robinson, 1967) to solve the normal equations (Robinson, 
1964). An error in EUREKA's algorithm was found 
(Montalbetti, personal communication); an updated 
version of EUREKA was used by the main program. 

In the study of short period reverberation removal, 
the trial trace was considered as a section of a time 


series composed of white noise, Cor Cyr Cor cee S307 
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Figure 4.1. Short period reverberation removal. 
Upper sleft: ethe strial torace. wUppersgrights ithe 
autocorrelation plot of the trial trace. 

Middle left: deconvolved trial trace - a=l, m=13. 
Middle right: deconvolved trial trace - a=3, m=ll. 
Lower left: deconvolved trial trace - a=4, m=10. 
Lower right: deconvolved trial trace - a=5, m=9. 
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Figure. ¢2 smsnOrtmpe! OdeLevetpe latvoneremo vals 


Upper left: the trial trace. 
autocorrelation plot of the deconvolved 


a=4, m=10. 
Middle left: 
Middle right: 
Lower left: 
Lower right: 


deconvolved trial trace 
deconvolved trial trace 
deconvolved trial trace 
deconvolved trial trace 


Upper right: the 


Ernale-trace 


m=9. 

m=l11l. 
m=10. 
m=10. 
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Figure 4.3. Long period reverberation removal. 
Uppers lGlUc stile tidal trace, Upper right: the 
autocorrelation plot of the trial trace. 
Miradd@le lerc: “deconvolved trial trace = G=16, m=28. 
Middlesrignts thesautocorle lation plot of the 
deconvolved trial trace - a=16, m=28. 
Lower left: deconvolved trial trace - a=16, m=27. 
Lower right: deconvolved trial trace - a=4, m=40. 
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convolved with a characteristic wavelet, the trial 
wavelet. Since the time series is the convolution of 
a minimum-delay wavelet with white noise, equation 
4.10 shows that the autocorrelation of the time series 
is the autocorrelation of the characteristic wavelet. 
Wien =the autocorrelation of the trial trace was taken, 
the multiple on the autocorrelation plot was disregarded. 
To avoid errors due to the multiple, the section of the 
autLocorrelatvom plot Wsedsgwasmrestrictedsto, lags, .0 to 
St 

Lhe Sewo Vartablesswenesoamand Min the following 
discussion, they will be written as a pair a-m; for 
example 4-10 would correspond to a=4 and m=10. Accord- 
ing to the criteria stated previously, the optimum case 
should’ be 6-8; that is, a=6, m=8, and atmtl= 15; a=6 
is the lag corresponding to the second zero crossing on 
the autocorrelation plot, while atm+tl= 15 corresponds to 
the length of the characteristic wavelet and to the 
Signi ta canteporeionecnmtne autocorrelation plot. The 
trials were taken by varying one of a or m while holding 
ehevother Constantsmeother trials held atm+l constant or 
else: he lorgeaam iiconstantye ther case corresponding to, zero 
delay least-squares inverse filtering. The best combina- 
tion was 4-10; a= 4 corresponds to the first negative 


maximum on the autocorrelation plot; the width of the 
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output ‘spike#is®also 4.° This’ case is plotted on the 
lower left of Figure 4.1, and its autocorrelation is 
prottedgon thewupper*rignteor Figure 472. This auto-— 
correlation plot shows the degree to which the auto- 
correlation approaches zero from lags Geto Orwlags 
Ae combar 

When atm+l was held constant at 15, increasing 
a caused the spike to become two sided and also added 
NOUS Reece recase = )-9 Hls plotted On the lower right of 
Figure 4.1. The supposedly optimum case 6-8 showed the 
spike as a two sided pulse with the second leg of 
greater amplitude than the first, and considerable 
noise was also evident. Decreasing a added noise; 
the cases 1-13 and 3-11 are plotted on the middle left 
and middle right of Figure 4.1, respectively. 

When a was held constant at 4, little change was 
noted when m was increased or decreased; the cases 4-9 
and 4-11 are plotted on the middle left and middle right 
OLtmrgurel4 278 respectively. 

When m was held constant at 10, increasing a 
caused the spike to become two sided and also added 
noise; decreasing a added noise. The cases 3-10 and 
5-10 are plotted on the lower left and lower right of 


Figure 4.2, respectively. 
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When a was held constant at 1, the best output 
was-ak-sotmt] =.15;,) the, cases: for m.increasing,or 
decreasing showed noise buildup, until noise was of 
the same amplitude as the signal. The case 1-13 is 
plotted.on the middle.left of Figure 4.1;. these cases 
correspond to zero delay least-squares inverse filter- 
eTiClie 

In the study of long period reverberation removal, 
the trial trace was considered as a section of a time 
series containing a reflection and a multiple of the 
Same amplitude and shape but opposite sign. In this 
case, the multiple on the autocorrelation plot was 
regarded as real but undesirable. 

The two variables q and m were unrestricted. 
According to the criteria stated previously, the optimum 
easeyshouldsbe 16-28: that is, a=16, m=28, and 
gintl= 45; ¢= lo isthe lag vat.whichsthesfirst points,ot 
theamull tiple. is seecnsgon thew~autocorrelationsplov, while 
atm= 44 is the lag at which the last point of the multi- 
ple is seen on the autocorrelation, and atmtl= 45 is 
themiength of ithe trial trace. sincetrialsawere, taken 
by varying one of a or m while holding the other 
constants wother trials held gtms) constant. wihe. best 
combination was 16-28; this case is plotted on the 


Middleslert ot Migure 4.3, and ts aucocorrelacion 1s 


apeso geod ii.t stlipest 26 eae 
ayodfii sersvai versipe-sSessl ySISb orSs oF DOgees209 
fayoms: noijevedxevey bolted paol 30 Ybucte) ema ah 
anit 6 tO HOijo=e s 2s bsrobienoo asw sosad Inttd oda 
sit Yo sigitium 6 Hae hoistoe lier 6 painisshos asixse 
-pie edizogqa Jud sqede Sas shb3itqes ames 
esw dolg noLtslorroses0s Sa7 no elgisium odd ,92680 
-Sitiegteebm J2¢ [ase 2s -bebseyem 
bsisiitesinu sisw m bas » esidelasv ows ath | ©) 
nuimitqo odd Vievotvexg heseoe prratiy) sia o+ potiraesA 


aiand ni 


brs: .8S/=o ~aL=n ,2t' ten :8S+al ed bisoda Sess | 


to jaiog text eft dotdw de psl oft ai 41 =p got = Lhe 
alidw ,Jolg heisalezxosctus edd so ese <i sigiv ium ema) 
~Edivm sit 20 tnieq tasI oi dotriw 35 pst oft 2: bb =e 


fe ae, be apa sxonodus add no aso at stg 


asin sale atkins oar -o0nr9 feixd edo to dépnsl othe 
eae ene 


50 


Rlottedson athe .middlesrighteot, Migure.4 236 asThis 
autocorrelation plot shows the degree to which the 
autocorrelation approaches zero from lags a to atm 
or lags 16 to 44. This case shows that the multiple 
has been reduced in amplitude by about one-half, but 
it also shows that a second multiple has been added. 

When atm+l was held constant at 45, increasing 
or decreasing a added low levels of noise which 
increased as a was further increased or decreased. 

The exception was the case 17-27 which is plotted on 
the lower left of Figure 4.3 and which shows a con- 
siderably higher level of noise. 

When a was held constant at 15 or 16, little 
change was noticed as m was increased or decreased 
except the case 16-27 which shows a considerably higher 
Level OL ino1se:. 

Whenem= was Held constant at 27, 26, 39, or 30, 
little change was noted; however, when m was 27, a 
considerably higher level of noise was present. 

One trial was done for the case 4-40; this case 
is plotted on the lower right ‘of Figure 4.3; 0=4 1s 
the same as in the optimum short period case. This 
trial showed that short period and long period reverber-— 


ations could be removed simultaneously. 
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4.4 Single Channel Predictive Deconvolution of 


Seismic Data 


The method of single channel predictive deconvolu- 
tion was tested on some of the seismic reflection data 
obtained at the University of Alberta in 1971. The 
data was chosen from the same shot on September 2, 1971 
as was used in Chapter Three for instrument deconvolution; 
five ihannels were filtered iby a five sto thirty! Hertz 
zero phase shift eight pole Butterworth bandpass filter; 
the time interval chosen was eighty-five samples long, 
from 94.125) £O\/./5 seconds7ss#the data fis plottedyon the 
left-hand side of Figures 4.4 and 4.5. 

The single channel predictive deconvolution main 
program using subroutine EUREKA was used; the program 
Operates on one channel at a time although its arrays 
are in multichannel form. The positive lag half of the 
autocorrelation is plotted in the center of Figures4.4 
and 4.5; a good degree of similarity between the auto- 
correlation plots may be seen. The first negative 
maximum on each autocorrelation plot occurs at lag 5, 
while the second zero crossing occurs between lags /7 
and 8. 

Since the short period reverberation removal 
tests indicated that the prediction distance a was the 


important parameter to determine, whereas the parameter 
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Figure 4.4. Single channel predictive deconvolution. 
Left-hand side: five channels of data from the first 


shot on September 2, 1971. Center: 
plots for the five channels of data. 


autocorrelation 
Right-hand side: 


five channels of deconvolved data with prediction 


distance a = 5 and atmtl = 14, 
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Figure 4.5. Single channel predictive deconvolution. 
Left-hand side: five channels of data from the first 
shotyons Septemberjizye 1071. Center:autocorrélation 
plots for the five channels of data. Right-hand side: 
five channels of deconvolved data with prediction 
distance a = 1 and atmtl = 14. 
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m had little effect and the parameter at+tm+l was of 
only slightly more importance, a set of tests was 

run holding atmt+l constant at 18 while varying a from 
imeCtomres. With otmtl = 18, two cycles of the auto— 
correlation plot were spanned. It was found that the 
changes between tests as a increased were very subtle; 
when a was 1 to 5, the main event was reduced to one 
leg, and when a was greater than 5, two or more legs 
seemed to be present. Thus, a was chosen to be 5, 
although it might have been chosen to be 4, 6, or 7 with 
Pucite. Signi ti cantachange. 

The next set of tests was run holding a constant 
At Oeewhd le .Setbbinguddmt le cual stony.) 1 jal 4, ph esac 
20. Little difference was seen between the cases for 
an equal to ll, 14, or 17. According to the 
previously jstated oxriterion yi ithe, autoconreiations 
are highly regular, atmrlis set to Span a plus one 
complete cycle. One cycle on the autocorrelation plot 
was 9 lags, so at+tmtl was chosen to be 14. This case is 
plotted on the right-hand side of Figure 4.4. 

A final test was, nune withsos sets at. 1 sand o@tmtl 
set at, 14. ‘This case, corresponds to wero delay least- 
squares inverse filtering and is plotted on the right- 
hand ‘side of Fagure- 4.5% “lt may be seen that the 


deconvolved data is somewhat sharper for this case than 
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for the case with a set at 5 and a+tm+tl set at 14. In 
both cases, the output data was filtered by a five to 
CoUGevenerezeczerTORDNise simi tecignt pole Butterworth 


bandpass filter. 


4.5 Multichannel Predictive Deconvolution 


Multichannel seismic reflection data may be 
deconvolved by means of a multichannel least-squares 
inverse filter which is determined by requiring that 
the sum of the squared errors between the desired outputs 
and the actual outputs for each channel be minimized. 
This requirement leads to a set of simultaneous linear 
equations, called the normal equations, each element of 
which is itself a square matrix (Wiggins and Robinson, 

Lo GS) 
Let the multichannel input be represented by X, 


where t is an integer; then, 


T 
Ss (Kip rXogrX gure Xz) ’ 4.36 


where T denotes the matrix transpose. The filter is 


a pie 


i -/7fh where each i rnin oy ete te ty uly ee er 


Oren ihiie is 
gxk matrix. The actual multichannel output is 


ak 
es (YoerYourYaur- ++ You) i 4.39 


The relationship between the input and the actual output 
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is given by 


ae f x,t fy t-1 + £oX-_2+ oie Pee ttm . 4.40 
The) désared multichannel output is 

on = tz Z; Z Zz Ae 4.41 

io eZ es 3 Cee ee Ct : ; 
Therefore, it is required that 

4 2 
ie ze: Elly,, — 2,,)7] aeAe 
i=l 


be a minimum. When the multichannel input represents a 
multivariate random stationary process, the normal 


equations may be stated in matrix form as 
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and 


in Onder to use equation 4243 to do multichannel 
predictive deconvolution, the number of input channels 
Muse equal Sthnesnumber ofsoutput, Channels, that @1s7k—10. 
Then, the set of cross-correlation coefficients may be 


specified as 
Gelr=oir The, i > =O IZ OLAS ny 4.46 


where q is the prediction distance. The set of filter 
coefficients fort, f5,----f, thus represents the multi- 
channel prediction operator. The multichannel prediction 


error is defined to be 
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4.6 Multichannel Predictive Deconvolution of Seismic Data 


The method of multichannel predictive deconvolution 
was tested on some of the seismic reflection data obtained 
atutheatniversityn ofeAlbertas inhie7 lease Thes-data was 
chosen from the same shot on September 2, 1971 as was 


used in Chapter Three for instrument deconvolution; 
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however, the four channels chosen are not the same as 
the five channels used to test the method of single 
channel predictive deconvolution. The data was filtered 
by ative to thirty Hertz zero phase shift eight pole 
Butterworth bandpass filter; the time interval chosen 
Woacmc COM mi =etOmo. OmSeCCONGG, stnie Gdatamis plotted ain 
Figure 4.6. The four channels were aligned by picking 
the initial trough on the event at 7.5 seconds; the 
shift from trace to trace was three samples. 

A main program was written which used subroutine 
AUGURY, which made use of subroutines OMEN, FORM, MAINE, 
MOVE, ZERO, BRAINY, and HEAT (Robinson, 1967). AUGURY 
does multichannel predictive deconvolution for a 
prediction distance a equal to 1; MAINE calculates 
the inverse of a matrix; BRAINY does multichannel 
convolution; and HEAT calculates multichannel auto- 
correlations and cross-correlations. 

With the alignment as specified, the data was 
deconvolved by multichannel prediction; it was found 
that three of the output traces were picked with peaks 
aligned but that the second trace from the top had 
troughs where the others had peaks, and vice-versa. 
One-half cycle of the dominant energy in the data was 
four samples; when the second trace was shifted by a 


positive or a negative lag of four samples, the traces 
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Figure 4.6. Four channels of data from the Peas eet tO 
on September 2, 1971. Each channel was filtered by a 
five to thirty Hertz Butterworth bandpass filter. 
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were all picked with positive peaks. The case for 
which the second trace was shifted by a negative lag 
of four samples was chosen for the remaining tests, 
and isis this lccase which vseplottedein Figure 4.6. 

A set of=tests wasexun holding g constant at 1 
and. Secuinggutmt | CQuaile COmsn. 4) oieewscn, LO: thus’ m 
ranged from 1 to 8. The case qtmt+l = 3 was very similar 
to the original data; the similarity between output and 
input decreased as atm+l increased. The cases qtm+l =9 
and 10 showed strong similarities between traces but 
few strong events and little similarity to the original 
data. The cases atm+tl = 6, 7, and 8 seemed to be a 
compromise; they showed some similarity to the original 
data, good similarities between traces, and strong events 
running across the traces. The case atmtl = 7 was judged 
to be the best; this case is plotted in Figure 4.7. 

As a comparison of the methods of single channel 
and of multichannel predictive deconvolution, the same 
four channels plotted in Figure 4.6 were deconvolved 
using the single channel main program which used sub- 
routine EUREKA. The prediction distance q was set equal 
to 1 and atmtl was set equal to 7; the output is plotted 
in Figure 4.8. The output from each method was filtered 
by a five to thirty Hertz zero phase shift eight pole 
Butterworth bandpass filter. It may be seen that the 


multichannel output in Figure 4.7 has slightly stronger 
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Figure 4).7/. 
Four channels of data deconvolved by multichannel 


prediction with prediction ‘distance @ = 1 ‘and 
OIL = 7. 


Multichannel predictive deconvolution. 
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Figure 4.8. Single channel predictive deconvolution. 
Four channels of data deconvolved by single channel 


prediction with prediction distance a = 1 and 
atm+l = 7. 
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events running through it than has the single channel 
output in Figure 4.8, but the similarity between the 

two plots is quite striking. This. may. suggest, that 

the method of single channel predictive deconvolution 

is as effective as the method of multichannel predictive 
deconvolution in defining major reflections. If this 

is correct, then the method of single channel predictive 
deconvolution is the preferable one to use, as it is 
considerably less expensive in terms of computer costs. 
In one application of multichannel least-squares inverse 
filtering of common depth point seismic reflection data, 
it was concluded.that the results were_no better than 
those obtained using single channel least-squares 
inverse filtering (Davies and Mercado, 1968); that 
conclusion appears to be in agreement with the conclusion 


reached here. 


besssiioxvatia aan 
svitpibsrgq fophedbisivm ze pinay C ere . 
eit 31 -enobipa ites a10fan patatiap at noboutowcns®: 
evidoibsig Lsnnsis sfpnia to horstem sit nord \498%100 aL 
at di 26 , Seu oF Sho slds19i9%g bas et TOstuLovacTsh. 7 
.eteon t9sugme> to ented aL avianeqee eeol yidsrebiencs 7 
setevai agveupa~teasel Llonnsdoisinum to noisesilggs sno al 
,69ab aoisoslisx otmetee toitog dsgsb nommco To patvos ize a 
aerie sej36d on stew etivass sft tans bebulougs aaw a | 
estsupe-—taangl Lenco olpaia pntes boatesdo seodt as 


jadt 4 (80€L ,ObsoxeM bac 2oivsd) pnitosli2? eezeval 
notewionos sit diiw tasmesips at ed oF pie eames Ae : 


BIBLIOGRAPHY 


Aipasian, J wl 968. . Spectral Behavior of Short Period 
Body Waves and the Synthesis of Crustal Structure 
in western Canada, M.Sc. thesis, Edmonton: 
University of Alberta, Department of Physics. 

Allsopp, Dik., Burke, M.D... and Cumming, G-Lisgy, 1972. 

A Digital Seismic Recording System, Bulletin of 
the Seismological Society of America, Vol. 62, 
pp. 1641-1647. 

Chandra .NON and Cumming _~G.b., 1972. Seismic 
Refraction Studies in Western Canada, Canadian 
Jourmia lant sarin ssciences » Vol a9 27pp. sl009=Lipoe 

Clowes, R.M., 1966. Deep Crustal Seismic Reflections at 
Near-Vertical Incidence, M.Sc. thesis, Edmonton: 
University of Alberta, Department of Physics. 

Clowes, R.M., 1969. Seismic Reflection Investigations 
at Crustal .Structure in southern Alberta, Ph.D. 
thesis, Edmonton: University of Alberta, 
Department of Physics. 

Clowes, .R.M., and Kanasewich, E.R., 1970. Seismic 
Attenuation and the Nature of Reflecting Horizons 
within the Crust, Journal of Geophysical Research, 
Wool eit, Bs OU ober 

Clowes, R.M., Kanasewich, E.R., and Cumming, G.L., 1968. 
Deep Crustal Seismic Reflections at Near-Vertical 


Incidence, Geophysics, Vol. 33, pp. 441-451. 


64 


s20beyHt Bo saduiaiae agar a9 iL 


AVOL \.dVD ,emimmed bas \.0.M yoxmud . 2.0 vqgoet ia i 


jo didelivd ,meseye pnibseves oimeie® Istbpid A 
,809 .foV \sottomA to yt9i502 iaytpofomaie® sat 
4 


-VhOL-I8I Jaq 

oimete2 .STe€l ,al.D ,paimmw> bas ,.4.M ,etbasdd 
tstbsaso ,sbsns9 aisie#eW ni esibuse nolitoasiea = 

-ROL{-e@0l .qq.,.@ .foV ,asonsic® dirs to Ismivot | : 

ts encizoslisd aimais@ LstanzD qeaad .300f ,.M.A ,eowold 


:nosnomba ,etesi+ .08.M ,soashionrt Lsbid1sy—te0K 


-GQ.d9 ,sdisdfA oxedjuo? mt suutounse Istestd to 
HSrodIA to yjtexevinU :notnomba ,aieon? . 
-abbayild to tnomttsaqsd 
oimete® .OTeL ,.A.d ,doiweesnsx bos ,.M.a ,eewolD 7 
anon izoH Limam 20 suvts sd+ brs aotteuns33A 
D 20 Lemwot \Fent) oft aitebw 
; ,20Ta-£e92 “aa ,e tov 7 
w eae peoreaae Be oe .29WOlD 
ror, eee 


a) 


65 


Cunmincy Gls. ja Ganland .G. De eandeVoZort.. Keen ooo. 


Seismological Measurements in Southern Alberta, 
Final Report, Contract AF 19(604)-8470, Project 
Vela-Uniform, Air Force Cambridge Research 


Laboratories, Bedford, Massachussetts. 


Cumming. G.li.., and Kanasewich, 2b -R.,-1966., Crustal 


Davies, 


Doh, 


BLES; 


Fuchs, 


Structure in Western Canada, Final Report, 
Contract AF 19 (628)~-2835, Project Vela-Unifomn, 
Air Force Cambridge Research Laboratories, 
Bedford, Massachusetts. 

E.B. and Mercado ,-E.Jo, 1968.5) Multichannel 
Deconvolution Filtering of Field Recorded Seismic 


Datal sGeophysics ,. VOl Soy Upp aw 22 


Ge, sana Puchns) kaw, e106 potatisticals Evaluation 


of Deep Crustal Reflections in Germany, Geophysics, 
MO ios 2 DD ie o/s 

REM. cance bashan jew Obs Crustal Character 
istics from Short Period P Waves, Bulletin of the 
Seismological Society of America, Vol. 58, pp. 
1681-1700. 

K., 1969. On the Properties, of Deep. Crustal 
Reflectors, Zeitschrift fiir Geophysik, Vol. 35, 


Done oe L4oe 


Ganley, D.C., 1973. .A Seismic Reflection Crustal Model 


Near Edmonton, Alberta, M.Sc. thesis, Edmonton: 


University of Alberta, Department of Physics. 


tev Ma TO - 7 

S08L \.4 ,Tosov Bas ,.d.9 basizsd ,.d.d josimeesl " 7 

‘sitedih miediguoa Al s$nomowenel Isdipolomated 

sSopor9 , 0 bG-(b00) CL BA 2581909 dxoqea tenke | 
deunseest epbixdms’? ootot 22° (abtinU-siov 
addseeudoseesM (Biotbea@ ,esitossiedsd 

fefeuxd .3a0f \. 8.8 \dolwoesaed Bas ,.1-9 .pnimud 
,J10qea [snid \absensS mietesW oi sivsovise 
,mrortiat-sleV soefotd ,t€es-(8Sd) el TA jopx12n05 
.a9ixodsiodsd dorssaceA spbhixndms) so1r0t 212A 

-stteennoseasM ,bioiSessi ~— ; 

[snesdoisziom .80@I ,.0.a ,6bsoreM bas ,.4.4 ,esived 
simeis2 bebroseh biertt to pnitesliaq nortulovnacosd 
.SST-L1T .gG ,f€ .1ov ,avieydgoad ,s7s0 

moiseulsva ([soitefssie .Td0L ,.4 ,efout bas ,.9 ,andod 
,aoteydqosd ,ynamied at enokjoeiiea Lasewi> qeed to 
.Tee-f2e .aq ,S€ Lov 

-tsgosishiD {seau1D .88CL ,.W.S ,mpnesa bas , Mea ,eilid 

efit io nitellud ,eeveW I Boizsd trode mort eoisel — 

qq ,8€ .LoV ,soLzemA to yYXeiso® Inoipolometo2 
.OOTf-18aL 

indewx qoo00 to setdzeqord sft HO .edeL ,.a% , ecient 
.2€ JoW tieyiqosd 22 StiTHoesieS axososl tes 
-OBI-€£1 .aq 

{shoM ledeuxrD noljoolten oimeis2 A .éT@l ,.5.a ,yslasd 
‘yHOsHOMDT ,ekeedt 152.M ,BsxedLA ,nosnomba teen 
-2dtaydd to Inemdtsqed ,stredfA to yJisreviad 


66 


Nawtyelies., *ea., L969" ~eThe barthn 4) Crust and Upper 
Mantle, Geophysical Monograph Series: American 
Geophysical Union, Washington, D.C. 

Heacock7 @).Gs,/ed., 1971) “The Structure and Physical 
Properties “OL the Earths Crust, Geophysical 
Monograph Series: American Geophysical Union, 
Washington, D.C. 

Hron,> Pao.  eCLi Leriaer OY Selection .Opebnases in 
Synthetic Seismograms for Layered Media, 
Bulletin of the Seismological Society of America, 
Vol Gln pp. 165-779". 

Hron, F.) andskanasewich, E.R. , L97L. sSynthetic 
Seismograms for Deep Seismic Sounding Studies 
Using Asymptotic Ray Theory, Bulletin of the 
Seismological Society of America, Vol. 61, 
|p © Seeded og By S02 ae LOK 

Junger, As, 1951. * Deep Basemént Reflections in Big 
Horm County 9Montana, Geophysics;sVol. 16, 
pe =99=510% 

Kanasewich, E.R., 1973. Time Sequence Analysis in 
Geophysics, Edmonton: The University of Alberta 
Press. 

Kanasewich,E.R., and Cumming, G.L., 1965. Near-Vertical- 
Incidence Seismic Reflections from the 'Conrad' 
Discontinuity, Journal of Geophysical Research, 


Vol.” 70% pp. "3441-3446, 


| oo om . ‘ 
T: a - 
pa - —— - tag 
- 7 
’ 
7 ba a 
; ee ' A : i eye — ee oa i ~ alt 
° art oe iC BVO Soa Je) = 2 sa 
% ot ‘ Bes ys? ~~ r : ~ >is y 7 
P| : - 
7 AR ie + ei 
° TOS Bh Lie OW Ft age L. Pa Va ; 
2 Sa od “ a han . . ——— 7 aot 
; bs nis — : 
_ 


Isorevida bre Batifouste itt ; «siboseet ‘ 


NOL Wbe ..5.4 
inoieyitgend \cauadehareai oa Yo Bosuges 
‘goin Isoteyriqosa nsots ms vsebsee Muesvont 
3.4 .dospaiitesn 
at 2seedd Ye noidésise xt BixsdinD iver 1.7 , 00% 
.S5ibsM Deteysd 40% amar pome fe sksoddieve 
,SditemA to. Yietooe Issipolomeise sd3 to nisoliva = 
Leyvagar -4q .f8 .1ev 7 
ditedsaye ,IVel ,.An doiwoesne® Bos .-4 , tote 
esibudae paibnave® vimets2 geo xot amex poma to8 
eis to mitelisd ,ytceA? eR oi soIqmyed prea 
18 ley ysoitaméA Yo yisksoe levipolameiss 
. O08f-eent . qq 
pid ai anditsolieg jadmeesa good . Leet {LA \tepaut 
ot .fov \adfeyiigosd ered VWIAGOD sol! | 
" -bfe-e88 aq 
nt eieyiena sige bam ai .-A.4 \dolwsesnsa 


67 


iNOpGi ia ue se ULakeyaC. lu eand Hart, ee J.) Sus. jo 6ce 
The Crust and Upper Mantle of the Pacific Area, 
Geophysical Monograph Series: American Geophysical 
Union, Washington, D.C. 

McDonelya! JeyeeDoOba, PA Mi uis, whol. Sengpusi, UR sims 
VansNostrand, R.G.;. and White, U7E., 1958. 
Attenuation of Shear and Compressional Waves 
ie Plercesolale, sy Geophysics, s V0) oy poe 2 le 
439. 

Meisner, R., 1967. Exploring Deep Interfaces by Seismic 
Wide-angle Measurements, Geophysical Prospecting, 
Vou. os upp wooo O—o ly ¢ 

OBrien, PN.S..,1969.. Some Experiments Concerning the 
Primary Seismic Pulse, Geophysical Prospecting, 
ViO1 SMU ty) Pp Leo 

Peacocke WW. #.ancd Lreitel . So. L969 es Predictive 
Deconvolution: Theory and Practice, Geophysics, 
Von SYNE a) slase nal elles eke J 

Richa Loss C ane Walker, a0. 8e oOo CosULenenceOn 
the Thickness of the Earth's Crust in the Albertan 
Plains of Western Canada, Geophysics, Vol. 24, 
Boner aoa. 

Robinson, E.A., 1954... Predictive Decomposition of Time 
Series with Application to Seismic Exploration, 
Ph.D. thesis, Cambridge: Massachusetts Institute 
of Technology, Department of Mathematics. 


Reprinted in Geophysics, Vol. 32, pp. 418-484. 


2sVaW Riera he ‘bas vere 


-I&b iqg \£S lov ,eotayiqesd , slede Pe - on 
| al : 
oimerea yd asosinstnt gost) palsolaqua Toel ,-8 ,wadetoM 
, PALISsqeord “fepteyilaood \ednomonves3M slpie-sbhiw a 
-Vie-882 .gq ,aL stev 9 
oid prinisono> esnemitegxd ee 381 2.14 ynoiaa'o 7 
, paisosazord tno haey sane Diteios yinmisd 
-The-Lhe .gg Wi ,fov 
evisothbord ,Rd@f ,.2 ,lodiei? Gas dia .tsossed 
veokeydgoss ,soisoserd bas yitosdT isoLIvLovacesa ~ 
»  ROERSEL .gq 48 toy 
0 Snemexverat .e2eL \4L.d aetieW bas \.9.7 .abeedoba 
nstisdhA oi a su .Paae ried atts small 


et 
ohVs fe ere deco rst 


0 on vst - pire so 


68 


Robinson, E.A., 1964. Wavelet Composition of Time Series, 
and Recursive Decomposition of Stochastic 
Processes, in Econometric Model Building: 

Essays on the Causal Chain Approach, Wold, H.O.A., 
ed., Amsterdam: North-Holland Publishing Company. 

Robinson, E.A-; 1966. ‘Collection of Fortran Pl Programs 
for Filtering and Spectral Analysis of Single 
Channel Time Series, Geophysical Prospecting, 

VO ela, soUDpLeMmeiiie 

Robinson, E.A., 1967. Multichannel Time Series Analysis 
with Digital Computer Programs, San Francisco: 
Holden=-Day, Inc. 

Schriever, W., 1952. Reflection Seismograph Prospecting - 
How It Started, Geophysics, Vol. 17, pp. 936-942. 

Somerville, P.G., and Ellis, R.M., 1972. “P=Coda Evidence 
for a Layer of Anomalous Velocity in the Crust 
Beneath Leduc, Alberta, Canadian Journal of 
Barth Scrences, Vol. 9, pp. 845-856. 

Sprenke, K.F., 1972. An Application of the P—-Coda Spectral 
Ratio Method to Crustal Structure an Central 
Alberta, M.Sc. thesis, Edmonton: University of 
Alberta, Department of Physics. 

Steinhart, J.S., and Meyer, R.P., 1961. Explosion Studies 
of Continental Structure, Carnegie Institute of 
Washington Publication 622, Washington, District 


of Columbia. 


cone ishom aidenenoot Ab - ysteoooes 

a «fy OB sbiow ,dssemqga nisdd Laie ot io sqseod ps 

-yasgnod patna nds bos { foh-ds104 sins Sams be 7 aa 

emetperd TEI agxtrod to iteattes _paet 4.8 3 \noenidos 
shpnie to etayienA [srjosq@ Sas pirinesLi4 103 
,pniistoegectd fesieydgosdD «sired omit? LorrisdD 
-inemelagque .bf .foev 

2faylenA esizse omit [sansdotsiumM veel ,.A.8 aeeniaen 
:oD2ibosia ase ,semetpord astugqmod Istipid sti 


ol , ysd=nebioH ; 


— pridocsqeoi® jasipometea fottoafisA .S2@eL ..W -asveltdo2 ory 
j 


~SS8-dF© .qq WE -LoV ,20fayiqosS ,bsdsisde 2t won ae 

POR sboo-F .SVEL ,.M.A, peilla bok, .0.4 \silivremea iM 

jel Sat ci yatbatav auoismonsh 26 3S5y¥sa 6 tot ’ 

Yo Lsntuol mpibensD <s27Sdifk ice Atsensa a 
.028-2bS -49 +2 ad aaa faend 

isisosge BuOI-F ors. to colgnottqas ha sates, a Pe, % simone | © 

faisneo nt sxvdsurs2) Eageuyd ce) boddoit. often 


he yi SPENT: 8 tnesoneR ee Se oad 


69 


SLetiiiaGe,. 0 toe, aNCwOMcll, Blin redo 6 lob ome mia ren 
Beneath the Continents, Geophysical Monograph 
Series: American Geophysical Union, Washington, 
Dee 

Vogel 24.7, .c¢., J9/ le. Proceedings of the Colloguaum ion 
Deep Seismic Sounding in Northern Europe held 
in Uppsala,on December 1, and 2, 1969,.. The. Swedish 
Natural Science Research Council (NFR), Stockholm, 
IPOH AES 

Wiggins, 2.-A..,4 and, Robinson, EaA.,¢ 1965... Recursive 
Solution to the Multichannel Filtering Problem, 
Journal of Geophysical Research, Vol. 70, pp. 


EUs s floyd bros Me 


oan 


aviewoent ..2aer a sang aided 
»moildard prixzos LEY fea FTO Lo 
-qq .0T .LoV _dloxs | 


- 7 = 7 - i : >a - 

ar Oy - 

‘cae pie “la ~~ ae 

a [fr Sage # ng Rew 

6 «gia S « ano a . ey a an at 
= | oa 

ry - cats ~— 


APPENDIX 


TRANSFER FUNCTIONS OF THE SEISMIC RECORDING INSTRUMENTS 


AND SEISMOMETERS 


In order to obtain the s-plane representation of 
the seismic recording instruments and seismometers, it 
is necessary to find an expression for each of the 
instruments. In the following paragraphs, the s-plane 
representation of the seismometers is derived. 

Each seismometer consists of a frame in which a 
mass m is suspended vertically by a spring with spring 
constant k and the motion of the mass is damped by a 
force proportionalgto velocity with damping constant q. 
If y represents the upward displacement of the frame, 
and gq represents the upward displacement of the mass, 
let x be defined as the relative upward displacement of 


the frame and mass, 


mg = kl j fee 


Where GpuSecthesaccelerationsol oravity.andal «is the 
SxtcistOnwoLetiecmspring, 4 Ute nes sci siGMeler Lrame 


moves y and the mass moves q, the equation of motion is 
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k(1+y-q) + a(y-q) - mg = mg . pe 
Using equation A.2, equation A.3 becomes 

kive-1qg)) 450 (y—q) = mq". ae 
Now, using equation A.1, equation A.4 becomes 


kx + ax = m(¥ - &) , Pe 
(Oe 


kx + ox + mé =my . A.6 


If the right-hand side is set equal to zero, the 


characteristic equation is obtained, 


If there is no damping, a = 0, and the solution of 
equation A.7 is sinusoidal motion with a natural 


frequency 


Thus, equation A.6 may be written 
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3 
x(s) +5 x(s) +5 x(s) =o V(s) A.10 


mw 
n On . 


Now, define a normalized complex frequency 


then equation A.10 becomes 


oO Zhe 
UGS lg SP iw, yeas gS) setae Ys Oy) Mare A.12 


the transfer function in S of output displacement 
for an input delta function in displacement is defined 


as 


Ty (S) = Mo ncrry Se TE aE - Axis 


The output voltage lefpropontional tosthebvelocity 


and is given by 
e(x) = ECy : A.14 


taking the Laplace transform of each side of equation 


A.14 gives 
Bits) =f-Gs xX .(S ae, AZLS 
or, in terms of the normalized complex frequency S, 


E(S) = “Gu 5 Gl Cs) c A.16 
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The transfer function in S of output voltage 
for an input delta function in displacement is, 


therefore, 
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The transfer function in S of an n-pole low pass 


ies ljeewe cali 


and the transfer function in S of an n=-pole high pass 


fiVeer"is 


T.,(S) a ee, Aw19 


where ForPy Fore er Fy are constants which describe the 
fl ltere’"A-bandpass’ tilter has a transfer function in S 
Witenes @ the=product Ot= the "Craiister LUNCtION Ino OL 
atlow pass’ filter*and the transfer function in S of a 
RiGhepass Liter. 


The constants in the s-plane representations of 


the seismic recording instruments and seismometers may 
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be found by examining the specifications of the circuits 


and their components. Then, these nominal values may 


be adjusted slightly so that the amplitude and phase 
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spectra calculated from the specifications match the 
amplitude and phase spectra determined experimentally. 
The constants in the s-plane representations of the 
seismic recording instrumentsand seismometers were 
determined according to this procedure by Cumming. 

In the following paragraphs, the parameter S 
equals iw/27. The function subprogram which expresses 
the s-plane transfer function of the seismic recording 
instruments used at the University of Alberta in 1969 


follows: 


COMPLEX FUNCTION TRANS (S) 
COMPLEX S 
TRANS= 
WG(S 7510171120 4S/ 410) 10/01 b42*574700 at (S/4 7008 ie eo) 
C THE TRANSFORMER 
Pate 7 213057010487 2130) 
C THE GAIN CONTROL OF THE PREAMP -- NOMINAL IS .120 
Te S169) /(186+87 169) 
C THE SECOND GAIN CONTROL 
IEUEYLBLITYVNIT OF8 77810) 
1ans/7e82 9) /(1.048/. 329) 
TFRIG7E 333) / (1.0487 .333))**2 
CATHE, LOW? CUT.» FPITGETERS 
TFL OVS OPSY57.3)) **4 


C HIGH CUT FILTERS’ -- NOMINAL VALUE OF 52.1 HZ. INCREASED 
C BY 1L0% 

1*({57+600)/7 (130757. 600) 
CG THE OUTPUT? TOeF .Mi“ TAPE -- NOMINAL IS .455 Ha. 

RETURN 

END 


The function subprogram which expresses the s=plane 
transfer function of the seismic recording instruments 


usedsat' the University of Alberta in 1970 follows: 


. Bxrdos¢ Sesiiq 2 - 
Ne de 
eit to ere <2 en: * hanna 


- aims ie ee ata a lites nel . 

2 tetamsisq sit , erlqsrpsrsa palwok Lot aig ai 
esesergze foidw mé1porgdus noidonu? one sno kuk i efseps: 
pniibiooss oimeies sit to noktsaue adbihicts snalq-e Say 
29el ni saTS ee te yetexsvial ofa tf beev einomutta " oy 
ated is 7 i 

EP 
‘\ 


_ 
- 


(2)}2HART WOLTOKUY SaLISMOD % 
2 XSrqMoo $ 


iehiad Celenanial s eatehade eel RG te (OL \eto, ry \UOE.Ae) [a 
SaMAOWRME , GHT 

(ObT.\et0. 1) \(OEL. Ney et 

OSL, 21. JAWIMOW -- 4MASAd aH? 70 JOoaTNOD WhaD SAT 

(ROL. \e45 .1)\ (eal. \2) 1 ‘a 

AJOHTHOS Ulad Ghose auF Dp 

(Ite, \a48. I)\(L1€.\8) 8k 

(OSE, \240, £)\ (O88. \e SE 

S*A( (EEE, \e+0. Ty) \ (£££. \2))*. 

SABEILG PUD WOl . we 

ae | EV A\8+0. t)\0.2 i” G 

astAanonI 8H D.S¢ 40 aUdAV JAMIMOU 7 BAMTGES 7 izK; P<. 


(woa.Nexo. t) wes a a> 


Bi aan “eI DAVIEUOM - =~ GaaT OM, si OT 


> 


i 7 1? . e 


Ip 


COMPLEX FUNCTION TRANS (S) 
COMPLEX S 
TRANS= 
o/e 10) it. 0rs/ 10) ) i 0/1 0+1  2*S/4700.4+(S/747002) #22) 
C THE TRANSFORMER 
ASUS /ia 20) CL eO+S/ ol20) 
C THE GAIN CONTROL OF THE PREAMP 
NS / bOI) (a0 eS / 2169 } 
C THE SECOND GAIN CONTROL 
HES ero Gl OF G7 =. 3 lel") 
Lz (S/o FAICCL S0ES An320)) 
TG Ceo) (Le OFS/ 2335) *e2 
CAT HESLOW CUTLYELETERS 
ere O/a(ls Ors 46 5.0)) ) #4 
Celie HIGH (CUE IMEGTERS (==¢NOMINAL, VALUE vES -52 18 HZ. 
LAGS) e207 CL. OFS/, 260) 
CHETHE “OUTPUT STAGE 
RETURN 
END 


The function subprogram which expresses the s- 
plane transfer function of the seismic recording 
instruments and seismometers used at the University of 


Alberta an L97L follows: 


COMPLEX FUNCTION TRANS (S) 

COMPLEX Ss 

TRANS= 

1((S/.10)/(1.0+S/.10) )*1.0/(1.0+1.2*S/4700.+(S/4700. ) **2) 

C THE TRANSFORMER 
1*1.0/(1.0+8/53.1) 

C THE PREAMPLIFIER 
1* (S/-100)/ (1.0+S/.100) 

Cath D.C. .OLP SEL FILTER 
fete oy (046/53. 1) 

C THE GAIN CONTROL 
1*105./(105.+105.*2.1105*S/50.+45.* (2.1105*S/50.) **2 
1+10.* (2.1105*S/50.) **3+(2.1105*S/50.)) **4) 

C THE FOUR POLE BESSEL FILTER WITH A FUDGE FACTOR FOR 

C FREQUENCY 
1* (S/7.5)**2/(1.04+1.4*S/7.5+(S/7.5) **2) 

C THE SEISMOMETERS 

RETURN 
END 
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Liemcne fourth Vast line of the subprogram and. its 
corresponding comment card are omitted, the function 
subprogram which expresses the s-plane transfer 
function of the seismic recording instruments used 
ai the Universicy of Alberta in 197INMis obtained. 

Figure A.l shows the amplitude and phase spectra, 
determined from the transfer function, of the. seismic 
recording instruments used at the University of Alberta 
vie OF. einem gures: Al WALZ, Als, and A.47eirequency 
is on a logarithmic scale, ranging from one-hundredth 
to one hundred Hertz, while amplitude, or gain, is in 
decibels and phase is in degrees. Figures A.2, A.3, 
and A.4 show the amplitude and phase spectra, determined 
trom, the transfer functions, of the seismic recording 
instruments used at the University of Alberta in 1970, 
of the seismic recording instruments used at the 
Uneversity of Alberta in 1971, and of the seismic 
recording instruments and seismometers used at the 


Dreversity of Alberta fin 1971, respectively. 


ee. 
beey etaomuxdent paibiongs ¢ 
bodtaido ef £0el nt imnade be 
,sisesge sesdq bas sbytilqms silt ayors ey, A 7 7 
oimaise sit to ,soisonut setensat ony mont ania na ir 
sitediA to ytiersvial ede ts been etnamessaak enkbsopey rt 
yonsupex® (8.4 DAB ,EsA 4S. A f.4 spaupit ai 00k dd Sa 
ap 


ddhoubayd-sno moti paipas: ,sisoe oimdsixspol ; no 2k 


at 2b \nkag Yo ,sbisilqms slinw ,stied Sorbapd sao OF | cs, 
,f.A 48.8 eorlpli -2spipsO mi ei Segng bas cho AG 


y 
Basimreseb ,sijosge sesdg bas sbustiqns sat wot dA bas os 
ar 
piibtove: aimetse Sdd io ,enotsatut tslensss ‘ead a a 
,OUET i Pree i te yjiersvinU sis 75 boex etoomeztenk aan 
4 ; 


6i+ ts boeu esnemuztens poidisset opmaige ott to 3 


7 aac] 


. 
piméise' sit to has ,IVeL ai stusdlA to eect - ane 


aig 36 boat eratemome ise Sas asnomsseai 
Vlevitesgesx ,fVel ar 4 


1 


PHASE (DEGREES) 
-36.00 36.00 108. 00 160.00 


-108.00 


-01 10.0 100.0 # 


1.0 
FREQUENCY (HERTZ) 


-10.00 


-20.00 


nw 
= 
io 
ea) 
_— 
O 
ug 
2 


S 


GAIN 
-40.00 


1.0 
FREQUENCY (HERTZ) 


Figure A.l. Amplitude and phase spectra of the seismic 
recording instruments used at the University of Alberta 
ry ARG 9. 
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Figure A.2. Amplitude and phase spectra of the seismic 
recording instruments used at the University of Alberta 
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Figure A.3. Amplitude and phase spectra of the seismic 
recording instruments used at the University of Alberta 
Bt 97. 
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Figure A.4. Amplitude and phase spectra of the seismic 
recording instruments and seismometers used at the 
University of Alberta in 1971. 
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